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We present the multispin coding for the nonequlibrium reweighting method of the Monte 
Carlo simulation, that was developed by the present authors. As an illustration, we treat 
the driven diffusive lattice gas model. We use the multispin coding technique both for 
the spin update and for the calculation of the histogram of incremental weights, which is 
needed in the calculation of nonequlibrium reweighting. All the operations are executed 
by the bitwise logical commands. 
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1. Introduction 

Monte Carlo simulations are used as standard techniques to investigate statistical 
mechanical properties of many-body systems. The use of efficient algorithms is 
important, and much effort has been dev oted to the development of new Monte 

Car lo alg or it hms , such as cluster alg 

orithmsC™ 

and extended ensemble methods 



For systems with discrete symmetry, a refined treatment is possible for the coding 
of the computer program. For instance, in the simulation of the Ising model, only 
one bit is required for storing the information of a single spin; a computer word 
can store the information of several spi ns a t the same ti me . The multispin coding 
technique developed by Bahnot et al. ^1 and Michael ^2 is based on this fact, 
and the Ising spins of 64 systems are put in a single computer word when using 
the 64-bit machine. The Monte Carlo spin flip process is executed with the logical 
commands, and 64 Ising systems are updated simultaneously with a single random- 
number sequence. Accordingly, the computation time is reduced remarkably. For 64 
systems, one may assign either systems having different parameters, for example, 
the temperature, the external field, etc., or systems having the same parameters but 
with different random number sequences. The multispin codin g has b een successfully 
used for the study of the Ising model on the regular lattices ■ l -*^*^' and that on the 
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quasicrystal El The three-state random Potts model together with the block-spin 
transformation was also formulated using the multispin coding The multispin 
coding was applied not o nly to the single spin flip dynamics but also to the Kawasaki 
spin exchange dynamics ^^^1 

The histogram rewcighting El jg a powerful method to investigate the equilib- 
rium properties with simulations; only simulation at a single temperature is requi red 
to obtain information for a range of temperatures. Recently, the present authors 1201 
have extended the idea of reweighting t o the c ase of nonequilibrium systems based 
on the Sequential Importance Samphngl^E3, A sequence of the micro-states, or a 
path, is generated, and the average over a single path is performed to calculate ther- 
modynamic quantities for a standard Monte Carlo simulation. In nonequilibrium 
reweighting, on the contrary, many paths are first sampled with a trial distribu- 
tion, and thermodynamic quantities are calculated based on the relative probability 
between the trial distribution and the target distribution. The relative probability 
is called "weights" in literature, which we shall use hereafter. The nonequilibrium 
reweighting technique was used for the nonequilibrium relaxation of t he I sing model, 
and the dynamical properties of the phase transition are discussed Moreover, 
the present authors 1^ appli ed the nonequilibrium reweighting to the study of the 
nonequilibrium steady states 1^^^ t he driven diffusive lattice gas model proposed 
by Katz, Lcbowitz and Spohn (KLS)ESlwas employed as an example of the system 
showing the phase transition to the nonequilibrium steady state. 

In this paper, we present the multispin coding of the nonequlibrium reweighting. 
It greatly saves com putation time. We pick up the driven diffusive lattice gas model, 
the KLS modelEi, to illustrate the multispin coding, but the formulation is general. 
Of course, it can be used for the Ising model, which is simpler. 

We organize the rest of the paper as follows. In Sec. |2 we briefly explain the 
KLS model. In Sec.|21 a short description of the nonequilibrium reweighting method 
is given. Sec. 0] which is the main part of this paper, gives the illustration of the 
multispin coding. In Sec. El we give the result of the calculation. The final section 
is devoted to summary and discussions. 



2. Model 

We use the KLS model l^iil for an illustration of the nonequilibrium reweighting 
method. This model consists of a half filled lattice gas on a periodic x Ly square 
lattice. Its Hamiltonian is given by 

n = -4 n.jTiej' (1) 

where the sum of lattice sites {ij} is over nearest neighbors, and the variable at the 
lattice site riij = 1 if the site is filled and riij = if the site is empty. For the Monte 
Carlo update, particles are allowed to hop to an empty nearest neighbor site with 
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the Metropolis rate, 

TpM^'W) = min[l, exp(-/3(A7^ - e^))] (2) 

where cr and a' are the system configurations before and after the hop, ATi. repre- 
sents the change in energy due to the hop, i? is a constant driving force, e = —1,0 
or +1 depending on whether the hop is against, orthogonal or along the direction 
of the drive, and {3 = 1/T is the inverse temperature of the heat bath. The Ly 
direction is taken as the direction of the drive. 

The KLS model exhibits an order-disorder second order phase transition. The 
ordered phase consists of strips of high- and low-density domains in the directio n of 
the drive. In the final steady state, a single strip of high density domain is formed 1^3 
The density profile along the direction of the drive can be regarded as the order 
parameter, and it can be defined as 




3. Nonequilibrium Reweighting 

In a Monte Carlo simulation, the system configuration changes with time. In each 
Monte Carlo step, a randomly chosen particle is attempted to hop to a nearest 
neighbor site and the attempt is either accepted or rejected. Let cti be the initial 
system configuration and let CT2 be the system configuration of the second Monte 
Carlo step, and 0-3 be the system configuration of the third Monte Carlo step and so 
on. In this way, a Monte Carlo simulation can be represented by a series of system 
configurations and we define a "path" of the simulation as follows, 

Xt = (o-i, ■ ■ ■ ,<Jt-i,crt) (4) 

Many paths are sampled and the thermal average of a quantity Q may be calculated 
by averaging over these paths, {Q{t))f3^E = (V'^) Sj=i *3(2?t )i where the sum is 
over all paths indexed by j, /3 and E are the inverse temperature and drive used 
in the Monte Carlo sampling process. The objective of reweighting is to calculate 
the thermal averages of Q at another temperature /?' and drive E' . This can be 
achieved by using appropriate weights w, 

n n 

m))p'.E' = Y.''iQ{xi)/Y.wi (5) 

i=i j=i 
To calculate the weights, the following steps are implemented: 

(1) Assume that a path xl up to some time t is sampled from a simulation at (3 
and E. 

(2) In the next Monte Carlo step, choose a pair of neighboring sites at random. If 
one of the two sites is empty, perform a Kawasaki exchange between the two 
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Tabic 1. Mapping of AH -tE to a four- bit pattern M. 
The transition rates (Eq. J^J) s-l^o given. 
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sites with the rate Tp^E{(^''' Wl)- represents the trial system configuration 
after the move. 

(3) The trial move is either accepted or rejected. In each case, we define an incre- 
mental weight Sw-' , 

a) If the move is accepted, crl^i = cr'"' and 5w^ = Tp> ,e' i'^''' Wi) /'^I3,e{<7''' Wi)- 

b) If the move is rejected, a^^^ = and Sw^ = [1 — T;3'^_e'((t'-'|o-j)]/[1 — 

(4) The weights at i + 1 arc given by multiplication: 

= wf X Sw-' (6) 

For each path j G {1, • ■ ■ ,n} these steps are repeated until a predetermined maxi- 
mum Monte Carlo time is reached. 

4. Coding Techniques 

We shall illustrate multispin coding for reweighting with infinite drive; general- 
ization to finite drive is straightforward. We use bitwise operations to update our 
system configurations and to calculate the weights. On a 64-bit machine, 64 systems 
are simulated in one run. For an L^; x Ly lattice, we use L^Ly words to represent 
simultaneously 64 lattices with each bit in the word representing one lattice site. 
Our multispin implementation is done systematically using truth tables. Firstly, we 
construct the appropriate truth tables for various tasks, and then wc implement 
them using bitwise operations in the computer. 

4.1. Multispin coding of the Kawasaki exchange 

In the Kawasaki exchange, a pair of neighbor sites is chosen and if one site is 
empty and the other site is filled, an attempt is made for the particle in the filled 
site to hop to the empty site with a rate given by Eq. In the coding im- 
plementation, wc should note that possible values of ATi. — eE for infinite drive 
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MaxA X exp(-12/3) MaxA x exp(-4/3) 
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1....1 





pi 


1..1 


1...1 


0....0 












MaxA X exp(-8/3) MaxA 

Fig. 1. The entries of pO and pi are assigned according to the transition rates. 

are oo, 12, 8, 4, 0, — 4, — 8, — 12, and — oo. For example, when the attempt is made 
to move a particle against the drive, ATi. — eE = oo and if the move is along 
the drive, ATi, — eE ~ — oo. If the move is orthogonal to the drive, e = and 
AH = —12,-8,-4,0,4,8 or 12 depending on the local particle configurations of 
nearest neighbors to the chosen sites. We map the values of ATi — eE into a four-bit 
pattern M = {M3, M2, Mi, Mq), which is illustrated in Tabled there, we also give 
the transition rate (Eq. 10)) for each ATi — eE. In the multispin coding, we introduce 
a dummy variable p, which takes to 3 in the present case, to implement the accep- 
tance and rejection procedure. Two arrays pO and pi are allocated for representing 
the probability of appearance of p. The array variable A for pO[A] and pi [A] takes 
to MaxA — 1. For the size of array, MaxA, we choose a large enough number, for 
example 2^"^. The values of each bit in pO and pi are set according to the transition 
rates as shown in Fig.^ Bits in the Ath element of the arrays, when put together as 
(pl[A], pO[A]), form a binary representation of the numbers 3(11), 2(10), 1(01) and 
0(00). The arrays pO and pi are related to the transition probability distribution as 
follows. Let A be a random number between and MaxA — 1, then the probability 
of getting the array elements with the bit pattern (pi [A], pO[A]) > 2 is exp(— 8/3), 
for example. We build the arrays pO and pi bit-by-bit and shuffle each bit indepen- 
dently. Note that the sequence of shuffling for each nth bit in pO and pi must be 
the same for both arrays. 

Then, the actual Monte Carlo procedure in the multispin coding can be sum- 
marized as follows; 

(1) Pick a pair of neighboring sites and perform a Kawasaki exchange 
trial. 

(2) From ATi - eE, determine M = (A/3, A/a, Mi, A/o). 

(3) Let A be a random number generated uniformly between and MaxA — 1. 

(4) Assign an acceptance bit A = 1 (acceptance) if 

8 X A/3 + 4 X A'/2 + 2 X A/i + A/o + 2 X pl[A] + pO[A] > 4. (7) 

Otherwise assign A = (rejection). Assignment of A can be implemented using 
Table El 

(5) Finally, the Kawasaki particle exchange is performed using XOR{®) operations 
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Table 2. Truth tabic for acceptance decision. A is a variable defined such that A = I repre- 
sents acceptance and A = represents rejection. 
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of A with both neighboring lattice sites. 

© A and n^'j' ^ n^'j' © A (8) 

4.2. Multispin coding of nonequilibrium reweighting 

Implementation of multispin coding for reweighting is much easier than that of the 
Kawasaki particle exchange. Possible values of incremental weights Swi are 



(9) 



The weights can then be written as a product of incremental weights, 

Wl = (5wi)''i(*)((Sw2)'''^^*H<^W3)'*'^^*n'^W^4)''^^*H'^W5)'''^^*H'5w6)''«^*^ (10) 

where h-[(t) ■ ■ ■ hglt) are the number of hits on the incremental weights Swi ■ ■ ■ Swq 
during the course of simulation from time 1 to t. Note that 6wo is irrelevant in 
Eq. H10|) . Hence calculation of weights has been reduced to accumulating histograms. 
In our implementation, we store the hits in a large array dhiik) — or 1 for i = 
1, • • • ,6 and fc = 1, • • • , for an arbitrarily chosen ig. Weights and histograms are 
updated only once every to steps using Eq. (fTn|l . and hl{t) = h{{t—to)+J2*^=idhi{k). 
For calculating dhi{k), we use the four-bit pattern M and acceptance decision A 
defined earlier. Then for each k = 1, ■ ■ ■ ,to, assign dhi{k) according to the truth 
table given in Table 13 To calculate the summation X]l=i dhi{k) in the multispin 
coding, we use the same routine to compute the tot al m agnetization or energy. The 
bit-counting routine, BITCNT, by Ito and Kanada can be used. The essential 
point is that we do not have to calculate the weights by multiplication at each time, 
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Table 3. Truth table for the mapping from {M, A) 
to dhi, i = 1, ■ ■ ■ ,6. 
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but we only need to calculate the histograms hi(t), which are obtained by integer 
operations. 

5. Result 

We performed simulations with reweighting on the KLS model with infinite drive. In 
Fig. 121 as an example, we show the time variation of the order parameter, Eq. ||2Jl, 
for 64 X 32 lattice with infinite drive. Simulations were performed at T = 3.160 
and data reweighted to nearby temperatures, T = 3.150,3.155,3.165, and 3.170 
(from top to bottom). Averages were taken over 4.096 x 10^ samples. We should 
note that the direct calculations at temperatures different from T = 3.160 gave the 
consistent results as the reweighted ones within statistical errors, which shows the 
effectiveness of the reweighting. The detailed analysis of the phase transition based 
on the finite-size scaling was reported in a separate paper 

We here mention about the performance of our multispin coding. The computa- 
tion time of the calculation with multispin coding for 64 systems is 40 to 60 times 
faster than that for independent 64 simulations without multispin coding for typ- 
ical system size. This is a rough estimate because such a comparison depends on 
the optimization of program. The overhead of computation time to calculate the 
hist ogra m hi{t) is the same order as the spin flip process, if we use the BITCNT rou- 
tinel^once per N single spin flip steps. This is the same situation as the calculation 
of the thermal quantities. In this way, the calculation of weights by the histogram, 
Eq. H10|l . reduces the computation time. More important is that the accumulated 
errors will be much reduced for histogram calculations compared to the calculation 
with the multiplication of incremental weights each time. 

6. Summary and Discussions 

We have formulated the nonequilibrium reweighting method that is convenient for 
implementing multispin coding. In addition, we have shown how multispin coding 
can be implemented for both particle hopping and nonequilibrium reweighting. All 
the computations arc performed in terms of integer with logical commands. As a 
result, a large increase of efficiency has been achieved. We would like to remark that 
the nonequilibrium reweighting method is general and may be applied to various 
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Fig. 2. Order parameter with infinite drive on the 64 X 32 lattice. From top to bottom, values of 
T are 3.150,3.155,3.160,3.165, and 3.170. 

models with different Monte Carlo updates. In this paper, we have treated the case 
of two component system with the Metropolis update. Nonequilibrium rewcighting 
can be applied to multi component systems, for example, g-state Potts models, and 
also to the heat-bath updat e. A generalization of the formulation is necessary, which 
will be discussed elsewhere EHI. 
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